Electro-magnetostatic homogenization of bianisotropic metamaterials. 
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We apply the method of asymptotic homogenization to metamaterials with microscopically bian- 
isotropic inclusions to calculate a full set of constitutive parameters in the long wavelength limit. 
Two different implementations of electromagnetic asymptotic homogenization are presented. We 
test the homogenization procedure on two different metamaterial examples. Finally, the analytical 
solution for long wavelength homogenization of a one dimensional metamaterial with microscopically 
bi-isotropic inclusions is derived. 
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I. INTRODUCTION 

The foundational claim of the metamaterial research 
community is the ability to engineer artificial materials 
with deliberate electromagnetic responses. However, de- 
spite considerable developments in the metamaterial field 
over the last 10 plus years, the question of how to quan- 
tify the electromagnetic response of a metamaterial is still 
very much an open one. The electromagnetic response of 
any material is quantified by the so called constitutive 
parameters of the materials (permittivity, permeability, 
etc.) The determination of these constitutive parameters 
in the long wavelength limit is the subject of this paper. 

One of the earliest attempts to assign values to the con- 
stitutive parameters of a metamaterial was the method 
sometimes referred to as as S-parameter retrieval^. This 
method lacks generality in that it assumes that the meta- 
material in question is isotropic and has no electromag- 
neto coupling, though it is recognized that it can be 
used to characterize anisotropic materials with diago- 
nal constitutive tensors. It has since been generalized to 
accommodate bi-anisotropic metamaterials obeying cer- 
tain symmetry requirements^ as well as reciprocal bi- 
isotropic (isotropic chiral) metamaterials^. In addition, 
there have been many other modifications of the origi- 
nal S-parameter retrieval method which are less useful, 
either because they are trivial or because they use the 
retrieval method in ways that are inappropriate consid- 
ering the assumptions the original method is based on. 
In addition to the fact that S-parameter retrieval is only 
applicable to highly symmetric metamaterials, it should 
also be noted that this method assumes no spatial disper- 
sion in the material it is applied to. This simplification 
affects many other assumptions implicit in the retrieval 
method regarding the symmetry of the constitutive pa- 
rameters and the boundary conditions of the material, 
limiting the validity of the S-parameter retrieval method, 
as well as introducing apparently non-physical charac- 
teristics to the retrieved constitutive parameters due to 
the presence of spatial dispersion^. Despite these limi- 
tations, S-parameter retrieval is today the most widely 
used method for assigning values to the constitutive pa- 
rameters of metamaterials. 

In addition to S-parameter retrieval, there have been 
several attempts at developing a more general homoge- 
nization theory for metamaterials^—. There has been 



limited success with these methods but none can be con- 
sidered a general theory for metamaterial homogeniza- 
tion. One trait that all of these methods have in common 
is that none of them are related to the classical method 
of homogenization known as asymptotic homogenization. 

Asymptotic homogenization, sometimes referred to as 
classical homogenization, is a old and established method 
for homogenizing differential equations that is limited to 
the long wavelength limit. For an introduction to the 
asymptotic homogenization method see Refs.QJ-QJ The 
fact that asymptotic homogenization only works in the 
long wavelength limit is the main reason that it has been 
largely ignored by the metamaterial theory community. 
Many of the more interesting metamaterial phenomenon, 
negative index of refraction for example, involve reso- 
nances which difficult to characterize using asymptotic 
homogenization. However, one advantage of asymptotic 
homogenization is that it can characterize very asymmet- 
ric and therefore anisotropic materials, something many 
other homogenization procedures cannot do. There have 
been a number of attempts to apply this method to meta- 
materials^—. In addition there is a similar homoge- 
nization method that involves placing a metamaterial be- 
tween two metal plates and calculating the macroscopic 
permittivity by treating the metamaterial as a capaci- 
to r 18 ' 19 . This method is less general than asymptotic ho- 
mogenization in that it requires certain symmetries in the 
metamaterial to be performed correctly. However, within 
this limitation the method is mathematically equivalent 
to asymptotic homogenization. In this paper we apply 
the asymptotic theory of homogenization to metamateri- 
als with microscopically bianisotropic inclusions. Our use 
of the method is distinguished from many of the previ- 
ously mentioned references in that we work with potential 
fields, as opposed to the electric and magnetic fields. 



II. ASYMPTOTIC HOMOGENIZATION 

A. The electromagnetic four-potential 

The normal implementation of asymptotic homoge- 
nization method is typically applied to a potential field. 
One example of this is the homogenization of electro- 
static permittivity tensors using the electrostatic poten- 
tial. Other examples include the homogenization of a 
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stiffness tensor using the displacement vector as the po- 
tential field, the homogenization of thermal conductiv- 
ity using the temperature as the potential field, and the 
homogenization of a permeability tensor using density 
as the potential field. Our first implementation of the 
asymptotic homogenization will use the electromagnetic 
four-potential, consisting of the electrostatic scalar po- 
tential Ao and the electromagnetic vector potential A. 
Following the standard method of asymptotic expansion, 
we assume that the four-potential can be represented as 
an expansion in terms of the unitless small parameter a 



A?(x) 
A«(x) 



Ag(x,y)- 
A°(x,y) 



<*Aj(x,y) 



-a 2 A 2 (x,y) + - 
-a 2 A 2 (x,y) + 



(1) 

Here y = x/o is a vector coordinate describing the 
rapidly changing values of the four-potential inside the 
unit cell and, x is a vector coordinate describing the 
slowly changing values of the four-potential on a larger 
spatial scale. The potential fields are periodic with re- 
spect to y, with periodicity given by the lattice constants 
of the unit cell. There is no requirement that the unit 
cell be cubic or tetragonal, only that it is periodic. The 
electric field and magnetic flux density are defined as 



E= -VAff, 
B = VxA". 



(2) 



1= I -A -1 , 

rh = — A 1 ■ C, 
q = pT l . 



(5) 



The differential operator acting on the four-potential 
in Eq. ((2]) is simply a four dimensional exterior deriva- 
tive (with the time derivative equal to zero) applied to 
a four- vector. In the traditional asymptotic homogeniza- 
tion, the equation of motion is an divergence (an interior 
derivative) applied to the constitutive matrix times the 
gradient (an exterior derivative) of the scalar potential. 
Similarly, our equation of motion is the interior derivative 
of the electric displacement and magnetic fields, which in 
turn are simply the constitutive matrix K times the ex- 
terior derivative of the four-potential 



= V-D 



p ■ (— VAg ) + f ■ (V x A") 



(6) 



= VxH = Vx[m- (-VAg 1 ) + q • (V x A a )] . 



Again we have assumed that the fields are static. One 
last detail is our choice of gauge 



V- A c 



0. 



(7) 



Here and throughout the rest of this paper we use 
Heaviside-Lorentz units. In Eq. ([2]) we have assumed 
that the fields are static. These fields are related to the 
electric displacement and the magnetic field by the mi- 
croscopic constitutive relations 



P [ 
rh q 

k 



(3) 



Here the constitutive matrix K consists of 3 x 3 tensors 
(p and q) and pseudotensors (I and rh). The microscopic 
constitutive parameters also vary according to both the 
large and small scale variables K = K(x, y) and are peri- 
odic with respect to y. These constitutive relations differ 
from the more typical constitutive relations 



C A 



(4) 



In the static regime the Lorenz gauge coincides with the 
Coulomb gauge. This gauge condition must be satisfied 
in addition to the field equation. 

As per the standard asymptotic homogenization 
method, we insert the field expansion of Eq. (JXJ) into the 
field equation Eq. (|6]) and gauge condition Eq. ([7]), sep- 
arating the resulting equations according to the order of 
a. The field equation proportional to a~ 2 is 



V y x 



K 



V y A° 

Vy X A 



0. 



(8) 



Here the operator \7 y only acts on the y variable of the 
potential field. The oT 1 order gauge condition is 



0. 



(9) 



This gauge condition combined with Eq. ([8j and the re- 
quirement of periodicity with respect to y implies that 
the 0th order fields are constant with respect to y or 
A° = A0(x) and A° = A°(x). 

or 1 order field equation is 



the relationship between the different microscopic consti- 
tutive parameters of Eqs. ^ and Q being 



V y x 



K ■ 



V^Aq V^Aq 
x A + V„ x A 1 



0. (10) 
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Here we have taken advantage of the fact that Aq and A 
are independent of y. We solve this equation by relating 
the first order four-potential to the derivatives of the Oth 
order four-potential 



Aj(x,y) = J2 



»=i 

6 



V x A°(x) 
- x A°(x 

VxAg(x) 



aoi(y), 



(11) 



a,(y). 



Here the subscript i on the 6x1 vector indicates the 
i'th component of that vector. We have introduced two 
new fields, aoi(y) and a$(y), which are typically called 
correctors. We ensure the periodicity of Aj and A 1 by 
forcing aoi(y) and aj(y) to be periodic with respect to 
y. With this representation we can satisfy Eq. ([TUf by 
solving the so-called cell problem 



Vy 

X7yX 



K 



V„A V^AJ 
X7 X x A 1 



V v x A z 



V^x 



j- , VajA^J 

y x x A + V y x A 1 



Y7 A 1 

VyA Q 



(17) 



= 0. 



Integrating this equation over the unit cell with respect 
to the y variable causes the first term to vanish due to the 
periodic boundary conditions. Using the representation 
for Aj and A 1 defined in Eq. (JTTJl, the integrated Eq. (fTTj) 
reduces to 



V x x 



K ■ 



-V^Aq 



= 0, 



(18) 



\y x x A J 

where the macroscopic constitutive matrix K is given by 



\VyX 



K 



-V y aoj 
V a x a.j 



0. 



(12) 



Here is a 6 x 1 unit vector pointing in the ith direction 



such that e,; • &j = Sij . 



The cell problem must be solved along with the a 
order gauge condition. 



V y ■ A 1 + V x ■ A = 0. 



(13) 



This equation must be true in order to satisfy the original 
gauge condition, but since it involves two fields we have 
additional freedom in how we satisfy it. The simplest 
choice is to use the Coulomb gauge for both fields, giving 
us 



and 



V a • A 1 = 0, 



V x ■ A° = 0. 



(14) 



(15) 



Equation (| 15[) is a gauge condition for the macroscopic 
potential A . Combining Eqs. (jT4"]) and (TTTj) we find 



0, 



(16) 



which is a gauge condition for solving the cell problem 
in Eq. (fT2"j). Equations (TT21 and (|TB|) must be solved in 
the unit cell of the metamaterial with periodic boundary 
conditions. We note here that the corrector fields aoi and 
a.i have units of length. Despite this, they can intuitively 
be thought of as the components of an electromagnetic 
four-potential. 

Finally, the a order equation is 



(K).. = - 



afy ei-K- 





f-V„aoA 


e 3 4 


Vv x a J/. 



(19) 



Here J Q d 3 y indicates an integral over the unit cell and V 
is the volume of the unit cell. The inverse volume term 
insures that the units are correct and that the averaging 
procedure returns the correct result when averaging an 
already homogeneous material. 



B. Dual scalar-pseudoscalar homogenization 

As an alternative to using the electromagnetic four- 
potential to perform asymptotic homogenization, we can 
take advantage of the fact that we are operating in the 
static regime, allowing us to represent both the elec- 
tric and magnetic fields as gradients of scalar and pseu- 
doscalar fields respectively. In addition to the previously 
used electrostatic scalar potential Ao, we now introduce 
the magnetostatic pseudoscalar potential Co- As before, 
we expand our fields in orders of the unitless small pa- 
rameter a 



A£ (x) = A° (x, y) + aAj (x, y) + a 2 A„ (x, y) + 
C£(x) = Cg(x, y) + ttCj(x, y) + a 2 C;](x, y) + 



The electric and magnetic fields are defined as 



E= -VA° 



H= -VC£. 



(20) 



(21) 



Because the fields are static we can disregard any vec- 
tor (or pseudovector) potential fields. This method will 
produce the same static constitutive parameters as the 
method outlined in the previous section, but with lower 
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computational cost. Instead of solving for a single com- 
ponent scalar field and a three component vector field, 
while enforcing a gauge condition, as was done in the 
previous section, we are now solving for two single com- 
ponent scalar (pseudoscalar) fields without the need for 
a gauge condition. 

The field equations that must be satisfied involve the 
divergences of the electric displacement and magnetic 
flux fields. 



= V-D= V- e • (— VAq) + £ • (— VCq) 



(22) 



= V-B= V- 



C-(-VA«) + A-(-VC£) 



By putting the expanded fields Eq. (I2U1) in to the field 
equation Eq. (|22p . and separating terms with the same 
order of a, we find a system of equations that must be 
satisfied. The a~ 2 order equation is 



c ■ 



/"-V„Ag\ 



(23) 



where C is defined in Eq. Q . This equation implies that 
the scalar (pseudoscalar) fields are independent of y or 
A°=A°(x)andC8 = C8(x). 



The a order equation is 



C ■ 



_V A° — X7 A A 



- v r 1 J 



= 0. 



(24) 



Here we have used the fact the both 0th order fields are 
independent of y. As is usual with the asymptotic ho- 
mogenization method, we represent the first order fields 
in term of the derivatives of the 0th order fields times 
correctors 



The a field equation is 



C- 



VyAg VajAj 



Va;Cn 



c ■ 



- V x Aq 

-V x Cq 



V^Aq 

v y c 



(27) 



= 0. 



Integrating this equation over the unit cell with respect 
to the variable y causes the first term to vanish due to 
the periodicity of the fields. Using the representation 
provided in Eq. (|25p . the integrated Eq. (f2"T)) becomes 
the macroscopic field equation 



c 



-v t a; 



o. 



(28) 



where the macroscopic constitutive matrix C is given by 



1 

V 



d 3 y et-C- 



— Vy&Oj 

-V„c 



(29) 



We note that the cell problem Eq. (|26|) and the formula 
for the macroscopic constitutive parameters Eq. ([29)) arc 
equivalent to those in Ref. [l6l However, our derivation 
follows the standard asymptotic expansion method and 
we have identified Co as the magnetostatic pseudoscalar 
potential. Our homogenization method uses the poten- 
tial fields whereas Ref. [l6| attempts to homogenize the 
electric and magnetic fields, and introduces two scalar 
fields without physical justification. 

Finally, a third formulation of the electro- 
magnetostatic asymptotic homogenization procedure is 
possible using vector and pseudovector potentials. The 
electric displacement and magnetic flux density would 
be defined as 



D = -V x C, 
B = V x A. 



(30) 
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co(*,y)=£u:s§ 



-V B A8(x? 



(25) 



c 0l (y)- 



This allows us to satisfy Eq. ([23]) by solving the cell prob- 
lem 



Vy 



c- 



-V y a 
Vj,c . 



= 0, 



(26) 



Here A is the electromagnetic vector potential used in 
SeclHlA and C is a magnetoelectric pseudovector poten- 
tial, the dual field of the electromagnetic vector potential. 
This method should return the same homogenization re- 
sults as the two previously described methods, but would 
be computationally expensive. It would require solving 
for two different 3 component vector (pseudovector) fields 
while enforcing two different gauge conditions. Though 
it is hard to imagine a situation where this method would 
be advantageous, it is still interesting to note that it ex- 
ists as an option. 



remembering that the correctors aoi(y) and Coi(y) are 
periodic with respect to y. Again, we note that the cor- 
rector fields aoi and Co^ have units of length, but can be 
thought of as electrostatic and magnetostatic potentials 
respectively. 



III. ONE DIMENSINOAL LAYERED CHIRAL 
METAMATERIAL 

As a first example of this long wavelength homogeniza- 
tion procedure, we examine a one dimensional metama- 
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terial consisting of periodically layered materials, some 
of these layers being microscopically chiral. 



y-2 



( 0.35a t 




! 0.35a ( 








e i Si 




£ i Si 


Ci w. 




Ci Ml 









►y-i 



The subscripts 1 and 2 indicate that the constitutive pa- 
rameters correspond to layers 1 and 2 (see Fig. [T]) . Notice 
that the microscopic constitutive parameters we are us- 
ing are reciprocal (e = e T , fi = fi T and £ = — ( T where 
T indicates the transpose) though this assumption is not 
necessary for asymptotic homogenization. Also, for sim- 
plicity we have made the chirality of the second layer in- 
dependent of frequency, despite the fact that physically a 
lossless chirality should be proportional to uj in the long 
wavelength limit. 



FIG. 1: Unit cell of a layered metamaterial. The metamaterial 
is periodic in the y Y direction with periodic lattice constant a. 
Translations in the y 2 and y 3 directions leave the geometry 
unchanged. The unit cell consists of two different layers, each 
with the constitutive parameters shown in Eq. pip . 

Fig. [1] shows the unit cell of the layered metamaterial. 
The isotropic constitutive parameters of the two different 
layers are 



£l = 


i. 


£2 = 


■5. 


Si = 


0, 


6 = 


2.85i, 


Ci = 
mi - 


0, 

-4: 


c 2 = 

M2 - 


-2.85i, 
— h 



(31) 



In the appendix of this paper we use the asymptotic 
homogenization procedure outlined in Sec. Ill Al to analyt- 
ically solve for the macroscopic constitutive parameters 
of a one dimensional layered structure. The resulting 
macroscopic constitutive parameters are 



C = 



(e x e± 0\ 

e|| f|| 

e || £|| 

Ci fx± 

C|| M|| 

\0 C|| 



(e/(e/i-C0) 



<e/(e/i - eO) W(e/i - SO) (S/(e/x - CC))(e/(e/x - £0) ' 

(fe^iO) 

(e/(6„ - CC)) W(e/i - SO) - (S/(eM - S0XS/(^ - so> ' 

(C/(^-S0) 

(e/(e/x - SO) W(eM - SO) - <S/(^ - S0)(C/(e/i - SO) ' 

W(^-so) 

(e/(e„ - SO) W(eM - SO) - <S/(^ - SO) Wfo* " SO) ' 

I 



(S), 
(0, 



(32) 



Here (X) indicates the arithmetic mean of X over the 
unit cell 



(X) = - Pd yi X( yi ), 
a Jo 



(33) 



where a is the lattice constant of the unit cell. In addition 
to the crystal being one dimensional, the only other as- 



sumption made in the derivation is that the microscopic 
constitutive parameters are bi-isotropic. 



After applying the formulas for the macroscopic con- 
stitutive parameters to the layered structure in Fig.[T]we 
calculate 
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FIG. 2: Isofrequency contours for (a) left handed and (b) right handed elliptically polarized waves calculated with an eigenvalue 
simulation of the layered structure shown in Fig. [1] (solid lines) and from the dispersion relation for a homogeneous medium 
with the constitutive parameters of Eq. (|34[) (circles). The labels on the solid lines are the normalized frequency ioa/c. The 
uia/c — 1.6 lobes on the left and right hand sides of Fig. [2jb) are left handed modes in the second propagating band. 




/'J 



3.81, 
-i4.75, 

i4.75, 
10.5, 



Cii 

/'I! 



2.20, 
i0.855, 
-i0.855, 
1. 



(34) 



Here j_ indicates the direction perpendicular to the layer 
interfaces or y l5 and || indicates the directions parallel to 
the layer interfaces or y 2 and y 3 . 

Looking at the long wavelength constitutive parame- 
ters in Eq. (1341) . we see that /Ltj_ is non- unity, despite 
the fact that no microscopically magnetic materials were 
present in the unit cell. This is due to the presence of 
the chiral layer. It is possible to create a magnetic re- 
sponse by adding bi-isotropic materials to a metamate- 
rial, though this comes at the cost of making the macro- 
scopic material bi- anisotropic. Second, notice that £j_ has 
a sign opposite that of £||, the same being true for dif- 



ferent components of C as well. Finally, the macroscopic 
constitutive parameters are reciprocal, which is to be ex- 
pected due to the fact that the microscopic constitutive 
parameters are reciprocal as well. 

To test these constitutive parameters we compare the 
dispersion relation of waves in a homogeneous effective 
medium with the constitutive parameters in Eq. (|34p . 
to the dispersion relation of waves in the inhomogeneous 
layered medium. Fig.[3]plots an isofrequency contour— 
of the layered metamaterial calculated using a finite ele- 
ment eigenvalue simulation (solid lines) as well as isofre- 
quency contours of a homogeneous medium with the con- 
stitutive parameters shown in Eq. (|34p (circles). Fig. [2] 
plots two modes, which are left and right handed ellipti- 
cally polarized. The isofrequency contours for the homo- 
geneous effective medium were calculated by solving the 
Maxwell equations in u and k space as a 6 x 6 eigenvalue 
problem. The eigenvectors correspond to the electric and 
magnetic fields and were used to determine the handed- 
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ness of each mode. This determination was done using 
the convention that a left handed elliptically polarized 
mode has an electric field at a fixed point moving coun- 
terclockwise around the wave- vector in time as seen from 
the point of view of someone the wave- vector is pointing 
away from. 

We see from Fig. [5] that the two dispersion relations 
agree very well when wa/c, k\a and k^a are all small. 
For large frequencies or wavenumbers, the dispersion re- 
lations diverge from each other. This is to be expected. 
The asymptotic homogenization method is only intended 
to be used in the long wavelength limit. For large fre- 
quencies unit cell inclusions can become resonant and 
at large wavenumbers spatial dispersion can become sig- 
nificant, both causing the long wavelength constitutive 
parameters to fail. 



IV. THREE DIMENSIONAL UNIAXIAL 
BIANISOTROPIC METAMATERIAL 

As a second example of bianisotropic asymptotic ho- 
mogenization, we consider a three dimensional crystal 
pictured in Fig. [3] 




FIG. 3: Unit cell of a three dimensional bianisotropic meta- 
material. The crystal has a cubic lattice with lattice constant 
a. Centered in the unit cell is a sphere of radius 0.3a consist- 
ing of an isotropic Tellegen material with constitutive param- 
eters given in the text. The material outside the sphere is a 
uniaxial dielectric with permittitivities provided in the text. 

The unit cell consists of cube with lattice constant a with 
a sphere at the center with radius 0.3a. The sphere is a 
nonreciprocal isotropic Tellegen material with the consti- 
tutive parameters 



e= 1.78, 

C= 2, 



£= 2, 
/i= 1. 



(35) 



The material surrounding the sphere is a reciprocal uni- 
axial dielectric with the constitutive parameters 



e_L = 


3.43, 


e ll = 


1.38, 


U = 


0.335, 


Z\\ = 


0.267 


Cx = 


0.335, 


Cn = 


0.267 


M-L = 


0.936, 


Mil = 


0.910 



e = c 



4, 
0, 



1.5, 
1, 



(36) 



where _!_ indicates the y 1 direction and || indicates the y 2 
and y 3 directions. 

The crystal shown in Fig.[3]was homogenized with both 
of the asymptotic methods presented in Sec. |TT] using the 
commercial finite element software Comsol Multiphysics 
3.5a. Copies of these homogenization simulations are 
available upon request by contacting the author using 
the email address preceding the references. The Upon 
applying either of the bianisotropic asymptotic homog- 
enization procedures outlined earlier in this paper, the 
macroscopic constitutive parameters are found to be 



(37) 



We can see that although the geometry of the unit cell 
is very symmetric, with several reflection and rotational 
symmetries, the macroscopic constitutive parameters are 
uniaxial and bianisotropic. This is because the micro- 
scopic constitutive parameters break the geometric sym- 
metries of the unit cell. The uniaxial dielectric surround- 
ing the sphere breaks several of the rotational symmetries 
of the unit cell, and the sphere consisting of the Tellegen 
material breaks all of the reflection symmetries. 

We test these macroscopic constitutive parameters by 
using them to calculate the dispersion relation of the 
crystal in Fig. [3J Fig. 2] shows two isofrequency con- 
tour diagrams, one calculated with a finite clement eigen- 
value simulation of the crystal in Fig. [3] (solid lines), 
and one calculated with the long wavelength constitutive 
parameters calculated with the asymptotic homogeniza- 
tion theory. From Fig. [4j we can clearly see that the 
two isofrequency contours agree in the long wavelength 
limit, but diverge for larger frequencies and wavenum- 
bers. The dispersion relation calculated from the macro- 
scopic constitutive parameters was obtained by solving 
the Maxwell equations in uj and k space as a 6 x 6 eigen- 
value problem. The resulting eigenvectors represent the 
macroscopic electric and magnetic fields. These eigenvec- 
tors have the same phase for the different components of 
the electric and magnetic fields, indicating that the eigen- 
modes are linearly polarized. Despite this, the eigenmodc 
polarizations do not in general lie upon the principle axes 
of the crystal. 

As a final note, looking at Fig. [H one can see that the 
contours for the eigenmode on the left side of Fig 2] seem 
to continuously merge with the contours of the eigenmode 
on the right side of Fig. 2] across the kia/n = ±1 bound- 
aries. A similar relationship exists at the k^aj^ = ±1 
boundaries. If one takes the two isofrequency contours 
in Fig. IU and tiles them in a checkerboard fashion, one 
sees that the contours for one eigenmode continuously 
merge with the contours of the other eigenmode as they 
cross the boundaries of the Brillouin zone. It is not clear 
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FIG. 4: Isofrequency contours for the three dimensional crystal shown in Fig. [3] Both isofrequency contours (a) and (b) are for 
different linearly polarized eigenmodes. The isofrequency contours are calculated with a eigenvalue simulation of the crystal in 
Fig. [3] (solid lines) and from the dispersion relation for a homogeneous medium with the macroscopic constitutive parameters 
given in Eq. (|37[1 (circles). 



why this happens, but it only seems to occur with a non- 
reciprocal bianisotropic crystal. 



Appendix A: Analytic homogenization of layered 
structure 



V. CONCLUSION 

In conclusion, we have presented two different imple- 
mentations of the asymptotic homogenization method for 
electromagnetic metamaterials with bianisotropic inclu- 
sions, including one method that handles a gauge con- 
dition. Because they are asymptotic methods, they are 
only valid in the long wavelength limit. We have tested 
the homogenization methods by comparing dispersion re- 
lations calculated from the resulting constitutive param- 
eters with the true dispersion relations returned by eigen- 
value simulations and found good agreement in the long 
wavelength limit. It is hoped that in the future these 
homogenization methods, particularly the four-potential 
method, can be generalized beyond the long wavelength 
limit to describe temporal and spatial dispersion effects. 



Here we analytically solve for the long wavelength con- 
stitutive parameters of a one dimensional layered struc- 
ture. We do so using the four-potential method presented 
in Sec. Ill Al The layered structure will have constitutive 
parameters that only vary in the y\ direction. They are 
constant in the yi and y% directions. This simplifies the 
problem considerably, allowing us to obtain an analytic 
solution. We assume that the constitutive parameters 
are bi-isotropic, thus all four tensors (pseudotensors) are 
symmetric under a spatial rotation and thus can be rep- 
resented as scalars (pseudoscalars) 
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The first step is to solve the cell problem in Eq. (fT7|) . 
To do this we first use the symmetry of the unit cell to 
simplify the cell problem. Due to the translational sym- 
metry of the unit cell in the 1/2 and 2/3 directions, the 
microscopic fields only depend on yi, or ao,(y) = aoj(yi) 
and aj(y) = aj(yi). Any derivatives with respect to yi or 
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2/3 vanish in both the cell problem Eq. (TT2|) and the cell 
gauge condition Eq. (TTS)) , The cell gauge condition be- 
comes d&ii/dy\ = 0. Thus all d&u/dyi terms in Eq. (fT2|) 
vanish. This leaves us with three out of the original four 
equations 



d 

dyi 
d 

dyi 



dyi 



Si, 



15a 



m5 2 i 



dyi 



(A2) 



d 

dyi 



mS 3i + q ( 

dyi 



Equation (|A2[) contains three different equations that 
must be solved with six different values for the index i, 
yielding 18 different equations overall. These equations 
must be solved with periodic boundary conditions on the 
yi domain. It should be noted that these three differ- 
ential equations are decoupled from each other and can 
be solved separately. Second, all three equations are for 
fields depending on a single dimension (yi). Finally, all 
three fall into one of three types of ordinary differential 
equation problems 



TABLE I: Solutions to Eq. 

all solved with periodic boundary conditions on y\ . Here 
■0 is a scalar field that only depends on yi and Ci_4 
are functions of y\ representing the various constitutive 
parameters in Eq. (|A2|) . The solutions to these ordinary 
differential equations are easily found and are shown in 
Table Q] Here the averaging operation () is defined in 
Eq. (03]). 



Now that we have the solutions to the cell problem we 
can calculate the macroscopic constitutive parameters. 
Using the simplifications mentioned above due to the one 
dimensionality of the unit cell, Eq. ([TO)) becomes 



dyi e* • K ■ 



( da, 0j 
dyi 

$2j 
S 3 j 

644 
da 3j 
dyi 

0&2j 



Si 3 \ 
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5j 



V dyi 



+ 5, 



(A4) 
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dyi 



= 0, 
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dyi 



C 2 +C; 



a 



dyi 



(A3) 



~ C4 d± 
dyi [ dyi 





d&oi 
dyi 


d&3i 
dyi 


da,2i 
dyi 


i = 1 


p(l/p) 








i = 2 





1 (m/q) m 
q (V<7> q 





i = 3 








1 (m/q) m 
q <!/<?> q 


i = 4 


1 (l/p) I 
p (l/p) p 








i = 5 





q(i/q) 1 





i = 6 








q(i/q) 1 



which when evaluated with the solutions to the cell prob- 
lem in Table |H yields 
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(l/q)(m/q) 



(1/9) 



(l/p) ' 11 
lm\ (l/p) (m/p) 

Pi (l/p) 



(l/q) 
(1/9)' 
(m/q) 
(1/9) ' 

911 = 



(1/9) 



(A5) 

Converting the K constitutive parameters in Eq. (|A4ft 
back into the C constitutive parameters gives us Eq. (|32[) . 
It is easy to see that in the absence of the bi-isotropic 
parameters (£ = C = 0) the macroscopic constitutive 
parameters reduce to the standard expressions. 
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